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Abstract 

The distribution of singular values of the propagation operator in a random medium is investi- 
gated, in a backscattering configuration. Experiments are carried out with pulsed ultrasonic waves 
around 3 MHz, using an array of 64 programmable transducers placed in front of a random scat- 
tering medium. The impulse responses between each pair of transducers are measured and form 
the response matrix. The evolution of its singular values with time and frequency is computed 
by means of a short-time Fourier analysis. The mean distribution of singular values exhibits very 
different behaviours in the single and multiple scattering regimes. The results are compared with 
random matrix theory. Once the experimental matrix coefficients are renormalized, experimen- 
tal results and theoretical predictions are found to be in a very good agreement. Two kinds of 
random media have been investigated: a highly scattering medium in which multiple scattering 
predominates and a weakly scattering medium. In both cases, residual correlations that may exist 
between matrix elements are shown to be a key parameter. Finally, the possibility of detecting a 
target embedded in a random scattering medium based on the statistical properties of the strongest 
singular value is discussed. 

Keywords: Ultrasonic waves in random media, Random matrices, Singular values statistics, Multiple 
scattering, Single scattering, Target detection 
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I. INTRODUCTION 



Wave propagation in a multiple scattering environment has been an interdisciplinary 
subject of interest in a huge variety of domains ranging, e.g., from solid state physics to 
optics, electromagnetism or seismology since multiple scattering can occur with all kinds of 
waves, whether quantum or classical. Among all areas of mesoscopic wave physics, some 
(like acoustics, seismology, microwaves) have the experimental advantage to offer control- 
lable multi-element arrays of quasi-pointlike emitters/receivers. In such a case, the propa- 
gation between two arrays is best described by a matrix, termed the propagation operator 
K. At each frequency, its coefficients kij correspond to the complex response between array 
elements i and j. Despite their diversity all practical applications of wave physics (communi- 
cation, detection, imaging, characterization...), have one thing in common: all the available 
information is contained in the array response matrix K. Once K is known, the rest is only 
post-processing. Therefore, in a random scattering environment it is essential to study the 
statistical properties of K, and their relation to field correlations, weak localisation, single 
versus multiple scattering. 

Previous works have been performed in a transmission context, whether it be for com- 
munication purposes [l-4] or scattering problems jll 0]. In this paper, we will consider 
backscattering configurations : the same array of N independent elements is used to trans- 
mit and receive waves. In that case, K is a square matrix of dimension N x N, and it is 
symmetric if the medium is reciprocal. We are particularly interested in the singular value 
decomposition (SVD) of the propagation operator, which amounts to write K as the product 
of three matrices: K = UAVt. A is a diagonal matrix whose nonzero elements Aj are called 
the singular values of K. They are always real and positive, and arranged in a decreasing 
order (Ax > A2 > ... > A^r)- U and V are unitary matrices whose columns correspond to 
the singular vectors. 
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It is now well known that in the case of point-like scatterers in a homogeneous medium, 

3, |8j, as long as 



long as 



each scatterer is mainly associated to one non-zero singular value of K 
the number of scatterers is smaller than N and multiple scattering is neglected 
So a singular value decomposition of K (or equivalently a diagonalisation of the so-called 
time- reversal operator KK*) allows the selective detection of several targets, each being 
associated to a singular value of K. This is the core of a detection method named DORT 



(French acronym for Decomposition of the Time Reversal Operator) 
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DORT has shown 



its efficiency in detecting and separating the responses of severa 



weakly heterogeneous media {l 11] as well as in waveguides 
in non-destructive evaluation 15], underwater acoustics [l6| 
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and in radar applied to forest environments 



scatterers in homogeneous or 



jl4j . It has found applications 



171 ]. electromagnetism 18M21] 



22-124 1 . 



In this paper we will deal with random scattering media, consisting of a large number 
(>> N) of randomly distributed scatterers, showing possibly multiple scattering between 
them. We will also address the issue of detecting a stronger reflector hidden in a statistically 
homogenous scattering medium, based on the strongest singular value Ai. Since the propa- 
gation medium is considered as one realisation of a random process, some general results of 
random matrix theory (RMT) may be fruitfully applied. 

RMT has been widely used in physics, statistics and engineering. The domain of appli- 



cations are numerous, ranging from nuclear physics 
networks telecommunications |l| or financial analysis 



251 ] or chaotic systems 26j to neural 
RMT predicts general be- 



haviours of stochastic systems as, for instance, determining the Shannon capacity for MIMO 
communications in random media jJ], 4} or statistical properties of highly excited energy lev- 
els for heavy nuclei 25J. Another direct application of RMT is to separate the deterministic 



and random contributions in multivariate data analysis 29^ 31] . 

Here, the experimental configuration we consider uses a piezo-electric array, with a finite 
(N =64) numbers of elements sending wide-band ultrasonic waves around 3 MHz in an 
a priori unknown scattering medium. The main issues we address in this work are the 

n 

applicability of RMT |1| to this experimental context, and its interest to establish a detection 
criterion based on the statistical properties of Ai. Since K is random, the relevant observable 
is the probability distribution function p(A) of its singular values. Recent experiments 32] 
indicate that in a multiple scattering regime, the distribution of the singular values is in 
good agreement with a simple law derived from RMT, the so-called "quarter circle law" 
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Pqc(A) — V4 — X 2 /tt (0 < A < 2) [j], [33I . In theory, the quarter circle law applies to 
square matrices of infinite dimensions, containing independently and identically distributed 
elements, with zero-mean and variance 1/N. 

Yet, from a physical point of view some of these assumptions do not hold: in a backscat- 
tering configuration the matrix elements are neither independent nor identically distributed, 
for several reasons (among others, reciprocity and weak localisation). In particular, the field- 
field correlations that may exist between matrix elements is a key parameter {24, 29 1. More- 
over when single scattering dominates, the spectral behaviour of K is shown to be similar 
to that of a random Hankel matrix (i.e, whose elements are constant along each antidiago- 
nal), and not a "classical" random matrix. This is due to the persistence of a deterministic 
coherence between single-scattered signals along the antidiagonals of the matrix K despite 
randomness, which does not occur in the multiple-scattering regime 
practical situation, the probability distribution function p(A) may differ significantly from 
the simple quarter circle law. Finally, in the last part of this paper we present theoretical 
considerations about the detection of a target embedded in a random scattering medium. We 
introduce a detection criterion based on the statistical properties of the strongest singular 
value Ai. This approach can be generalized to target detection in noisy environments. 



321 ] . Therefore in a 



II. EXPERIMENTAL PROCEDURE 



The experiment takes place in a water tank. We use a N-element ultrasonic array (N = 
64) with a 3 MHz central frequency and a 2.5-3.5 MHz bandwidth; each array element is 
0.39 mm in size and the array pitch p is 0.417 mm. The sampling frequency is 20 MHz. The 
array is placed in front of the medium of interest. The first step of the experiment consists 
in measuring the inter-element matrix (see Fig. [T|). A 100-/is-long linear chirp is emitted 
from transducer i into the scattering sample immersed in water. The backscattered wave 
is recorded with the N transducers of the same array. The operation is repeated for the N 
emitting transducers. The response from transducer % to transducer j is correlated with the 
emitted chirp, which gives the impulse response hij(t). A typical impulse response is shown 
in Figj2j The N x N response matrix H(£) whose elements are the N 2 impulse responses 
hij(t) is obtained. Because of reciprocity, hij(t) = hji(t) and H(i) is symmetric. We take as 
the origin of time (t = 0) the instant when the source emits the incident wave. 
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FIG. 1: Experimental setup. A 64-element array is placed in front of a random medium at a 
distance a. The whole setup is immersed in a water tank. The inter-element response kij(T,f), 
around the time of flight T and at the frequency /, is measured. It contains contributions of single 
and multiple scattering paths whose lengths belong to the interval [R — Ar/2;R + Ar/2], where 
R = cT/2 and Ar = cAt/2. Examples of a single scattering path (labelled s, dashed line) and of 
a multiple-scattering path (labelled m, continuous line) is drawn. (X S ,Z S ) are the coordinates of 
the scatterer involved in path s. (Xm\ Zm^) and (Xm \ Zm) are the coordinates the first and last 
scatterers along path to. 

A scattering medium is essentially characterized by its scattering mean- free path l e , and 
its diffusion constant D. If the scattering path length within the medium is larger than 
l e , multiple scattering can predominate. This is expected to happen at late times in the 
scattered signals hij(t). 

The impulse response matrix H(t) is truncated into short time windows in order to keep 
the temporal resolution provided by acoustical measurements and study the transition from 
a single scattering to a multiple scattering regime. The time signals hij(t) are truncated into 
overlapping windows : %(T,t) = Kj (T - t)W R {t) with W R {t) = 1 for t G [-At/2 , At/2], 
Wfi(t) = elsewhere. The value of At is chosen so that signals associated with the same 
scattering path within the medium arrive in the same time window. The detailed calculation 
of At is given in Appendix [A] In our experiments, we have typically At ~ 30 periods. At 
each time T, the kij form a matrix K. A short-time Fourier analysis is achieved by a discrete 
Fourier transform (DFT) and gives the response matrices K(T, /) at time T and frequency 
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FIG. 2: Impulse response hn(t) measured in the forest of rods (see Sec lIII Al) . 



/. The numerical SVD of each matrix K is performed and yields N singular values Aj(T, /) 
at each time T and frequency /. 

The next step consists in studying the distribution of singular values. It should be noted 
that we only have access to one realisation of disorder: the scattering medium is fixed, there 
is no ensemble averaging. Experimentally, the ensemble average can only be estimated by 
an average over frequency, and/or an average over time. The universal results provided by 
RMT are based upon the assumption that the elements of the random matrix have zero 
mean and a variance of 1/N. The first condition is easily met, assuming that kij(T,f) is 
a superposition of scattering contributions with a phase uniformly distributed between —tt 
and +7r. In order to fulfill the second condition and compare experimental and theoretical 
results, K is renormalized into K with normalized singular values Aj 



Once this renormalization is achieved, the experimental singular values distribution can 
be investigated. First, we form the histogram H(X) of the whole set of renormalized singular 
values Aj(T, /), taken at every rank i, time T and frequency /. Bins of this histogram are 
the intervals [mw; (m + l)w], where w is the bin width and m is a non-negative integer. 

denotes the number of singular values \(T, /) contained in the same bin as A. An 
estimator of the probability density function of the singular values is obtained: 




(1) 



P(A) 



H{\) 



(2) 



nw 
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n is the total number of singular values (n = N xnrxii/, tit is the number of time windows 
considered and n/ the number of frequencies at which the DFT of fey is achieved). At early 
times (cT — 2a<l e , with c the wave speed in the surrounding medium) multiple scattering 
can be neglected, whereas at later times (cT — 2a » l e ) multiple scattering dominates. In 
the following of the study, the theoretical singular values distribution predicted by RMT 
will be confronted to the experimental estimator p, both in multiple and single scattering 
regimes. 



III. MULTIPLE-SCATTERING REGIME 



A. Experimental configuration 



Here, we study the multiple-scattering regime. To that aim, we use a random-scattering 
slab consisting of parallel steel rods (longitudinal sound velocity cl = 5.7 mm//js, transverse 
sound velocity ct = 3 mm/ps, radius 0.4 mm, density 7.85 kg/L) randomly distributed with 
a concentration n=12 rods/cm 2 . The frequency-averaged elastic mean-free path l e is 7.7±0.3 
mm between 2.5 and 3.5 MHz 34j. The distance a between the array and the scattering 
sample is 25 mm. The slab thickness L is 40 mm. Even in such strongly diffusive media, 
single scattering occurs for small times of flight. Here, we will consider times of flight T 
larger than 70 /is, corresponding to scattering path lengths more than 7 mean free paths. 
Under these conditions, single scattering can be neglected and K only contains multiple 
scattering contributions. The whole experimental procedure described in SecJTT] is achieved 
and a set of renormalized matrices K(T, /) is obtained. The time-window length At is 
directly deduced from the calculation shown in Appendix |A] and is set to 10 ps. 



B. Experimental distribution of singular values 

We perform singular value decompositions of K(T, /) and obtain the estimator p(A) as 
defined in Eqf2J According to RMT, if the matrix coefficients are complex random 
variables independently and identically distributed with zero mean and a variance of 1/N, 
then the asymptotic {i.e . for N — > oo) density function of its singular values is given by the 
quarter circle law I3| • As we can see in FigfSl the experimental distribution of singular 
values is very far from this prediction. The major reason for this discrepancy is that the 
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FIG. 3: Experimental distributions of singular values of K (black disks) and Kt (white disks) are 
compared to the quarter circle law (black continuous line). The error bars are ± two standard 
deviations. 

coefficients kij are not independent. The experimental matrix K exhibits strong correlations 
between neighbour entries. They are measured by the correlation coefficient T m 



where the symbol < . > denotes an average over the variables in the subscript, i.e time T, 
frequency / and source/receiver pairs The integer m = i — j represents the distance 

between sources or receivers, in units of p (the array pitch). FigjU(a) clearly points out 
a strong correlation between neighbour entries, with a coefficient I\ = r_i ~ 0.5. The 
physical origin of these correlations will be detailed in Sec JIII Dl We will also show how to 
incorporate these correlations into the theoretical model for p(A). 

In this particular case, correlations can be simply removed by considering only one in 
two elements. From the initial matrix K of dimension N x N, a truncated matrix K< 
of dimension N/2 x N/2 is built, considering only elements of odd or even index. K t 
no longer exhibits short range correlations at emission and reception. The experimental 
distribution of normalized singular values is now much closer to the quarter circle law (see 
Figj3]). Nevertheless, a slight disagreement remains between the theoretical and experimental 
curves, specially at the vinicity of A = and A = 2. 
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FIG. 4: (a) Normalized correlations T m (real part : continous line; imaginary part : 
dashed line) estimated from experimental results, (b) The experimental distribution 
of singular values (black disks) is compared to the theory taking into account cor- 
relations (grey continuous line) and to the result of the numerical simulation (white 
disks) described in Sec lIIIDl The error bars are ± two standard deviations. 

In the next paragraph, we show that in the multiple scattering regime the variance of the 
matrix entries is not the same for all pairs (i, j), contrary to the assumption of identically 
distributed random variables. The consequence of this deviation from RMT assumptions 
will be discussed. 

C. Variance of coefficients kij. 

The signals kij(T, f) at a time T and frequency / correspond to the sum of partial waves 
that reach the array in the time window [T — At/2;T + At/2]. They are associated with 
multiple scattering paths whose length belongs to the interval [R — Ar/2; R + Ar/2], with 
R = cT/2 and Ar = cAt/2. An example of such paths is drawn in FigfTJ The response 
kij(T, f) can be decomposed into a sum of partial waves associated with the N q paths. In 
a 2D configuration, under the paraxial approximation and assuming point-like transducers 
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and scatterers, kij(T, f) can be expressed as 



^•(T,/)oc^fi, 



cxp 



jk ( zf + Z, 



r(2) 



9=1 



7(1) 7(2) 



exp 



.1 ; 



X, 



(1) 



2Z 



cxp 



jk 



Xj 



X, 



(2) 



2Z 



(2) 



(4) 

where fc = 2iif fc is the wave number in the surrounding medium. The index g denotes the 
q th path which contributes to the signal received at time T. (^Xq 1 ^, Zq 1 ^ and Z^ 2 "* j 

are respectively the coordinates of the first and last scatterers along the path q. B q is the 
complex amplitude associated with path q, from the first scattering event at yXq 1 ^, zip j 
until the last one at (^Xq 2 \Zq 2 ^. From the central limit theorem, we expect that the 
coefficients fcy are gaussian complex random variables. 

The renormalization (EqJTJ) leads to a set of matrices K. If the k^ were identically dis- 
tributed, the renormalized coefficients kij would be complex random variables with variance 
of - = 1/N. Actually, Figj5]^a) shows that of - is not uniform over the pairs The 





(a) (b) 

FIG. 5: (a) Mean intensity Iij of coefficients kij, normalized by jr. Results have 
been averaged over all times T and frequencies /. (b) Iij as a function of i — j (an 
average is achieved over couples such that (i — j) =constant). 



intensity /j 



ki 



averaged over time T and frequency / is displayed for each 



source/receiver couple It serves as an estimate of the variance ct?-. The diagonal 

elements ka exhibit a doubled variance compared to the off-diagonal elements. Moreover, 
the variance of off-diagonal elements tends to decrease when the distance \i — j\p between 
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the source and the receiver increases. From these results, it appears that the variance of the 
matrix entries, instead of being the same for all pairs (z, j), could be modelled as 

2^ 



No* = N 



kij 



Sv + f(\i-j\) (5) 



where 5 is the Kronecker symbol and / a slowly decaying function such that /(0) = 1. 
The doubled variance of diagonal elements compared to the off-diagonal ones is well 



explained by the coherent backs cattering effect 
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- |37| . It corresponds to an enhancement (by 



a factor of 2) in the intensity of waves scattered in the backward direction. This phenomenon, 
also known as weak localization, originates from a constructive interference between a wave 
traveling along a multiple-scattering path and its reciprocal counterpart: it appears when 
multiple scattering occurs, as long as the reciprocity symmetry is preserved. When source 
and receiver are identical (i = j), any path q whose first and last scatterers are of coordinates 
(^Xq 1 ^, Zq \ and (^Xq 2 \ Zq 2 ^j, interfere constructively with its reciprocal counterpart whose 
first and last scatterers are located respectively at yXq 2 \Zq^J and ^Xq~\Zq\. As a 
consequence, on average, the diagonal elements ka have a variance twice that of off-diagonal 
elements. Actually, this effect could also affect the matrix elements close to the diagonal. 



The typical width of the coherent backscattering peak is k ^/ D ^ 2 =j^ 38l . l39| , which has to 
be compared to the array pitch. In this experimental configuration, we have T > 70 fis, 
D ~ 4mm 2 / '/is 40], a = 25 mm hence — -, a < 0.16 mm, whereas the array pitch 



ky/D(T-2a/c) 

is p = 0.417 mm. Therefore the enhancement due to coherent backscattering is strictly 
limited to the diagonal elements. Finally, upon renormalization (Eq. [I]), taking the coherent 
backscattering effect into account amounts to rewrite the variance of k^ as afj = (see 
Appendix IB"1) . 

Now, one can wonder if the doubled variance of diagonal elements ka has any influence 
on the distribution of the singular values. It is shown in Appendix [B] that if the matrix 
dimension N is large enough, the influence of the coherent backscattering peak can be 
neglected as long as it is limited to the diagonal elements (y k J DT < Pj , provided that the 
matrix is renormalized according to EqJT] We will therefore consider that the coherent 
backscattering enhancement has no significant effect on the distribution of singular values. 

The doubled variance along the diagonal of K is not the only deviation from the identical 
distribution assumption. Indeed, if we focus on the dependence of afj as a function of (i—j) 
(see Figl^b)), we observe that it is not uniform: it decreases with the distance x — \i — j\p 
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between the source % and the receiver j. This is due to the progressive growth of the diffusive 
halo (i.e the mean intensity) inside the medium. It is not instantaneous but depends on the 



diffusion contant D 



41] . In a near-field configuration (a << Np, Np being the array size), a. 



would decrease as exp ^— wt)- And as soon as y/DT becomes significantly larger than the 

array size (Np), ofj can be considered as constant 42j. In a far-field situation (a » Np), 

afj is always constant, as long as % 7^ j. The experimental case we study corresponds to an 

intermediate configuration between near and far field. 

A solution to deal with the variations of afj is to compensate for them, afj depends 

2 

only on time T and \i — j\ (EqJSJ). Thus we can estimate of. by averaging kij(T,f) over 
frequency / and elements which are separated by the same amount m — \i — j\ : 



a 2 (T,m) = ( %(T,/) ) (6) 

\ / f,{(i,j)\m=\i-j\} 

where a 2 (T,m) is the estimator of afj. Then, at each time T, matrix entries fcy are normal- 
ized once again as: 

Sg(T, /) = M^,/) (7) 
\l^(T,\i-j\) 

A set of compensated matrices K C (T, /) is built from K. One can show that these matrices 
satisfy the identical distribution property : the variance of elements is now constant over 
all pairs Once this operation is achieved, the new distribution of singular values can 

be investigated and is plotted in Figj6]^a). However, the two sets of matrices K t and 
lead to similar singular value spectra. It means that in our case the variations of afj with 
— are not sufficiently large to significantly modify the distribution of the singular values. 
Nevertheless, in other experimental configurations, the non-uniformity of afj could have a 
stronger influence. For instance, in a near-field configuration, the growth of the diffusive 
halo is more visible and the distribution of singular values would exhibit stronger deviations 
from the quarter circle law, specially at early times. 



D. Influence of correlations on the singular spectrum of K 

As we can see in FigJD^a), the experimental matrix K(T, /) exhibits correlations between 
adjacent entries. There are two reasons for that. First, there is a mechanical coupling 
between neighboring array elements. Second, the wave recorded on the array can be seen 
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(a) (b) 

FIG. 6: (a)Influence of the non-uniform variance of kij'. The experimental distribu- 
tions of singular values for Kt (black disks) and (white disks) are compared to 
the quarter circle law (black continuous line), (b) Influence of the symmetry of K: 
the result of the numerical simulation (white disks) is compared to the quarter-circle 
law (black continuous line). The error bars are ± two standard deviations. 



as the radiation of a spatially incoherent source with width W (the diffuse halo inside the 
multiple-scattering medium) observed at a distance a. The Van Cittert-Zernike theorem 43|, 



44J states that the typical coherence length of the wave field is \a/W (in other words, 



the waves radiated by a finite-size incoherent source see their coherence length increases as 
they propagate). These two effects result in a short-range correlation between the scattered 
signals recorded on the array. In the experimental situation we investigated, the residual 
correlations are limited in range to adjacent elements, both in emission and reception, as it is 
shown in Fig3^a). They are mainly due here to the mechanical coupling between neighboring 
array elements, the coherence length of the diffuse wave-field becoming rapidly smaller than 
the array pitch in our experimental configuration. 



Sengupta and Mitra[29| have investigated theoretically the influence of correlations on 
the distribution of singular values in the case of large- dimension random matrices. Their 
model assumes first that (kij^j = 0. The correlation between two coefficients ku and kj m is 
expressed as 

(k u k*m) = N- 1 ^^ (8) 
where the symbol < . > denotes an ensemble average. C and D are N x N matrices. We will 



13 



refer to them as the correlation matrices. Based on EqJHl Sengupta and Mitra predict the 
distribution of singular values using diagrammatic and saddle point integration techniques. 
Only the eigenvalues of C and D are required to obtain the singular values distribution of 
K. 

We now apply the Sengupta and Mitra's approach to our experimental configuration. 
Here, the coefficients c^- and dij are identical because of spatial reciprocity : C = D. We 
estimate Cy by the correlation coefficient T^j defined in Eqj3] : 

Cij = Fi-j (9) 

The coefficients c^- form a matrix C which is an estimator of the correlation matrix C. C is 
a Toeplitz matrix : its coefficients only depend on % — j. The application of Sengupta and 
Mitra's approach to our experimental situation embodies two assumptions : 

• The correlation between two signals received (emitted) by the same element / and 
emitted (received) by two elements i and j does not depend on the element I but only 
on the distance \i — j\. 

• The correlation between two signals emitted and received by two pairs of elements (i, I) 
and (j,m) is given by the product of correlations at emission {T^j) and at reception 

(r«-m)- 

The first assumption requires the random medium to be statistically invariant by translation. 
The second one is made for analytical tractability and is verified in our case. 

Once the correlation coefficient T m is measured experimentally at emission/reception 
(see FigfJJa)), the correlation matrix C is estimated by C (EqJH]). Eigenvalues of C are 
calculated numerically and incorporated in Sengupta and Mitra's model in order to obtain 
an estimation of the singular value distribution for K. 

In FigH](b), the theoretical result provided by Sengupta and Mitra's method is compared 
to the experimental singular value distribution but also to the result provided by a numerical 
simulation. It consists in generating numerically a matrix P whose elements are circularly 
symmetric complex gaussian random variables with zero mean. Then, a matrix Q is built 
from P, such that 

Q = C^PC^ (10) 
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One can show that the matrix Q exhibits the same correlation properties at emission and 
reception as the experimental matrix K. A histogram of singular values is obtained by 
achieving the SVD of Q, renormalizing its singular values according to EqJTJand averaging 
the result over 2000 realizations. The agreement between the numerical and theoretical 
distributions of singular values is perfect, which illustrates the validity of the theoretical 



method given in [29|. Taking the correlations into account improves significantly the agree- 
ment between theory and experiment. 

As was mentioned earlier, when analysing the experimental data we crudely suppressed 
the correlations between neighbours by removing from the matrix most of its elements. One 
can object that it is silly to throw away the data, instead of trying to make the most of 
it. It is true that the correlations could have been cancelled out by applying the following 
operation 

R = C 1/2 KC 1/2 . (11) 

R has uncorrelated elements, while keeping the same size as K. This procedure amounts to 
"whitening" the data. The singular value distribution of R would then fit the quarter-circle 
law. However, this method would be relevant only if the correlations were exclusively due to 
the experimental device (here, a mechanical coupling between neighbouring transducers). As 
mentioned previously, there also exists a physical correlation due to the intrinsic coherence 
of the diffuse wave-field. And in Sec JIVI we will wee that there may be another form of 
correlation, a long-range correlation, that is typical of single scattering. In the general case 
of an unknown medium, single and multiple scattering coexist, and all correlations (short- 
and long-range) must not be cancelled out. Indeed, applying EqQj] blindly would affect 
equally single and multiple scattering contributions, whereas their statistical properties differ 
(see Sec ilVI) . This is particularly important for target detection in a random medium (see 
SecjV]), where one tries to detect the single-scattered echo from a target drowned in multiple 
scattering, based on the statistics of the singular values. 

In addition to short-range correlations between neighbouring elements, the matrix K is 
symmetric because of spatial reciprocity. This property results once again in a deviation 
from the quarter circle law theorem which assumes independent matrix coefficients jl(. We 
will consider here the case of the truncated matrix K f as defined in Sec lHI Bl in order 
to study the effect of symmetry independently from the presence of correlations between 
adjacent entries. Another numerical simulation is performed. It consists in generating a 
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symmetric matrix P of dimension y x y wnose elements are circularly symmetric complex 
gaussian random variables with zero mean. A histogram of singular values is obtained by 
computing the SVD of P, renormalizing its singular values according to EqJTJand averaging 
the result over 2000 realizations. The distribution of singular values provided by numerical 
calculations does not exactly match the quarter circle law (see FigJBfb)). More precisely, 
the deviation appears mainly in the vicinity of A = 0. The reason for this deviation is shown 
in Appendix By and large, the symmetry of K induces additional correlations between 
diagonal elements of the autocorrelation matrix KK^, which do not exist when the random 
matrix is not symmetric. These residual correlations are shown to be of the order of 4p. In 
our experimental configuration (N/2 = 32), they are quite low but sufficient to induce a 
slight deviation from the quarter circle law especially near A = which can also be observed 
in the experimental results (FigJBk). 

As to the other end of the singular value spectrum, an interesting theoretical result 
is that the distribution of singular values should remain bounded (for N » 1) even in 
presence of correlations. If correlations are ignored, the maximum value predicted by the 
quarter-circle law (which would apply if the matrix elements were independent and iV >> 1) 
is Xmax = 2. Here, taking T m into account, theoretical results show that the singular 
values cannot exceed A max = 2.5 (see FigH^b)). However, note that a slight difference 
between experimental and theoretical results remains near \ max (see Figs]3j H£b) &HD^a)). 
A possible origin for this small discrepancy is the presence of recurrent paths, i.e., multiple- 
scattering paths for which the first and last scatterers are identical or located in the same 
resolution cell. A partial wave associated with a recurrent path would display the same 
statistical behavior as single scattered waves, similarly to what happens in the coherent 
)ackscattering phenomenon, where recurrent paths tend to diminish the enhancement factor 



451 ] . Therefore, as the single scattering contribution results in a singular value distribution 



which is not of bounded support (see Sec JIVp . recurrent scattering may be responsible for 
the slight disagreement between experimental and theoretical results near \ max . However, 
the remaining disagreement between theory and experiment is small, and it is difficult to go 
beyond speculations about its origin. 



1(3 



E. Conclusion 



It appears that in the absence of spatial correlations between matrix entries kij, the 
experimental distribution of singular values is close to the quarter circle law. Nevertheless, 
some deviations from the assumptions generally made in RMT have been pointed out. First, 
the elements of K are not identically distributed. The reason for that is the inhomogeneous 
distribution of backscattered intensity, due to the coherent backscattering enhancement as 
well as the decreasing of the diffuse halo with the distance between the transmitter and the 
receiver. Another cause of deviation has been pointed out numerically: the symmetry of the 
array response matrix, which induces additional correlations between matrix coefficients. 
All these effects become negligible as long as N » 1. Yet they slow down the convergence 
of p(\) to the quarter-circle law. 



In 

Mitra 



;he presence of correlations between entries, the method proposed by Sengupta and 



291 ] allows one to calculate the distribution of singular values. As a preliminary con- 



clusion, in the case where multiple scattering dominates, the normalized response matrix 
K falls into the general scope of RMT and a theoretical prediction of the singular value 
spectrum can be achieved. Particularly, an upper bound A max for the singular values can be 
calculated. As it will be shown in the next section, this is no longer true if single scattering 
dominates. Note also that this agreement with classical results of RMT only holds in a 
diffusive regime. Indeed, if we approach Anderson localization (kl e ~ 1), one can expect 



that interference effects, such as hot spots 
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471 ] or loops 
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should strongly affect 



the distribution of singular values. Yet in our experimental configuration, we are far from 
strong localization (kl e ~ 100). 



IV. SINGLE-SCATTERING REGIME 



In this section, we consider the case of a weakly scattering medium. The sample under 
investigation is now a slab of gel (composed of 5% of gelatine and 3% of agar-agar), with 
thickness L ~ 100 mm and a mean free path l e ~ 1000 mm. In such conditions, the 
multiple scattering contribution is negligible. The array-sample distance a is 60 mm. The 
experimental procedure and the numerical processing are performed as described in SecjII] A 
typical example of a matrix K(T, /) measured experimentally is given in FigJT^a). Contrary 
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to the multiple scattering regime, the array response matrix K exhibits a deterministic 
coherence along its antidiagonals, although the scatterers distribution in the agar gel is 
random. We briefly recall the origins of this phenomenon 32| and its impact on the singular 




J m 
(a) (b) 

FIG. 7: Experimental results in a single-scattering sample (agar gel), (a) Real part 
of the matrix K(T, /) obtained at time T=265 (is and frequency /=3.1 MHz. (b) 
Real(white disks) and imaginary (white squares) parts of (3 m are shown as a function 
of m and are compared to the real(continous line) and imaginary (dashed line) parts 
of Eqd 



values distribution. 

As before, the signals kij(T, f) at a time T and frequency / correspond to the sum of 
partial waves that reach the array in the time window [T — At/2;T + At/2] except that 
only single-scattering paths are now considered. The isochronous volume is defined as the 
ensemble of points that contribute to the backscattered signal at a given time. It is formed 
by a superposition of ellipses whose foci are transmitter i and receiver j. In a far- field 
configuration, we approximate the isochronous volume by a slab of thickness Ar = cAt, 
located at a distance R = cT from the array and parallel to it. For simplicity but without 
loss of generality, we also assume that the reflectors as well as the array elements are point- 
like. 

In a 2D configuration, under the paraxial approximation, kij(T, f) can be expressed as 



hj{T, }) OC 2^ A* ex P 



d=l 



jk 



[xj - X d 
2R 



exp 



jk 



(Xj - X d ) 
2R 



(12) 
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The index d denotes the d th path which contributes to the signal received at time T. Xd 
is the transversal position of the reflector associated with this path and the amplitude Ad 
depends on the reflectivity of the scatterer. 

Let us express ky as a function of variables (x» — Xj) and (xj + Xj): 



kij{T, f) oc exp 



jk 



K X{ Xj 



^ A d exp 



jk 



[xj + xj - 2X d 
4R 



(13) 



deterministic term random term 

The term before the sum in EqfTjJ] does not depend on the scatterers distribution. On 
the contrary, the term on the right does; hence it is random. As a consequence, single 
scattering manifests itself as a particular coherence along the antidiagonals of the matrix K, 
as illustrated in FigfTl Indeed, in a given sample, along each antidiagonal (i.e. for couples 
of transmitter i and receiver j such that i + j is constant), the random term of EqfTjJ] is 
the same. So, whatever the realization of disorder, there is a deterministic phase relation 
between coefficients of K located on the same antidiagonal. It can be written as: 

(mp) 2 



0n 



k 



i—m,]+m 



(TJ) 



ku(T,f) 



exp 



jk- 



R 



(14) 



Note that this essential result is valid independently of the scatterers configuration, under 
two conditions: single scattering and paraxial approximation. The parabolic phase depen- 
dence along each antidiagonal predicted by EqJT4]is compared in Figff^b) with the coefficient 
f3 m obtained experimentally at time T = 265 fis and frequency / = 3.1 MHz along the main 
antidiagonal {i.e. i = 32). Theoretical and experimental results are in a very good agreement. 

In order to investigate independently the effect of the deterministic coherence along an- 
tidiagonals of K, other correlations that may exist between lines and columns of the matrix 
K have to be removed. These correlations are measured by assessing the normalized corre- 
lation coefficient T m (Eqj3]). T m (not shown here) spreads until \m\ = 2 for this experiment. 
Thus, the initial set of matrices K(T, /) has been truncated by keeping only one in three 
elements. Note that only correlations between lines and columns of K are cleared. The 
deterministic coherence along antidiagonals is long-range and hence is not removed by this 
operation. As before, the truncated matrix K t (T, /) is renormalized (EqfTJ. The experimen- 
tal distribution of singular values is shown in FigJHJ clearly, it does not follow the quarter 
circle law, even though correlations between neighbouring entries have been removed. This 
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FIG. 8: The experimental distribution of singular values /5(A) in the single-scattering case (white 
disks) is compared to the quarter circle law pqcW (dashed curve) and to the the Hankel law 
Ph (A) (continuous curve) which has been calculated numerically. The error bars are ± two standard 
deviations. 



is due to the deterministic coherence of single scattered signals along the antidiagonals of 

To our knowledge, this kind of random matrix whose antidiagonal elements are linked 
with a deterministic phase relation has not been yet investigated theoretically. But its 
properties are close to those of a random Hankel matrix, whose spectral behaviour has been 
studied recently 49(]. A Hankel matrix is a N x N square matrix whose elements belonging 
to the same antidiagonal (i + j = constant) are equal. Let {a p } be a sequence of 2N — 1 
complex random variables identically and independently distributed with zero mean. The 
Hankel matrix R built such as = ai +J _i, is said random. We have checked numerically 
that a random matrix whose antidiagonal elements are linked with a deterministic phase 
relation displays the same singular value distribution as a Hankel random matrix, provided 
that all elements of this matrix are zero mean and have the same variance. In the literature, 
Bryc et al. 49( have proved, for normalized random Hankel matrices, the almost sure 
weak convergence of the distribution of eigenvalues to a universal, non random, symmetric 
distribution of unbounded support. We will assume that this convergence property also 
applies to the singular values distribution of a random Hankel matrix. In the following, 
Pi? (A) will denote the asymptotic singular values distribution of a random Hankel matrix 
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and will be referred to as the Hankel law. To our knowledge, no analytical expression of the 
Hankel law has been found yet and only a numerical calculation can provide an estimate of 
p#(A). In FigJHl the experimental distribution of singular values of K t has been compared 
to the Hankel law. PhW is provided by a numerical generation of random Hankel matrices. 
The agreement between both curves is satisfactory. An important feature of the Hankel 
law is its unbounded support. As a consequence, in the single scattering regime, the first 
singular value has no bound, contrary to the multiple scattering case for which the first 
singular value could never be higher than A max in the asymptotic limit (N — > oo). 

V. PROBABILITY DISTRIBUTION OF Ai 

So far, we have studied the distribution of singular values as a whole, in relation with 
the importance of single or multiple scattering within the medium. We now focus on the 
strongest singular value, Ai. Knowing its probability distribution function is crucial for 
applications to detection. Imagine a random scattering medium, in which a target (i.e., 
a stronger reflector) may be embedded. Once K is measured, if the strongest normalized 
singular value Ai is above a certain threshold a, we will conclude (with a certain probability 
of error) that there is indeed a target. In order to assess the performance of a detection 
scheme based on a SVD of the propagation matrix, we need a reliable model for the pdf of 

V 

Assume the occurrence of a target at a depth R. The array response matrix K(T, /) at 
the time of flight T = 2R/c can be written as : 

K(T,/) = K T (T,/) + K R (T,/) (15) 

K T (T, /) is the contribution of the direct echo reflected by the target, and K R (T, /) the 
response of the random medium, which may include single scattering, multiple scattering 
contributions as well as additive noise. K R (T, /) can be seen as a perturbation (not nec- 
essarily small) of K T (T, /). As usual, the matrix is renormalized according to EqJTJ For 
simplicity, we will assume that short-range correlations have been removed, so that the 
relevant probability density function of the singular values of K R should be the quarter- 
circle law (for N » 1 in the multiple scattering regime) and the Hankel law (in the single 
scattering regime). 
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In a first approximation, K T is of rank 1 (actua 



ly. 



or a resonant or a large target K T may 
have more than one significant singular value |50M54l|). K R is random, and its normalized 
singular values have a probability density p R (X). Depending on the scattering properties 
of the medium (particularly, its mean-free path), at time T = 2R/c, p R (X) may follow the 
Hankel law (if single scattering dominates), the quarter-circle law (if multiple scattering or 
additive noise dominate) or a combination of both in intermediate situations. Let Af denote 
the highest singular value of K R (T, /). As a detection threshold based on the first singular 
value of K is needed, the statistical behavior of Af has to be known. If the singular values of 
K R (T, /) were independent from each other, then the distribution function F R of Af would 
be simply given by the N th power of the distribution function F R (X) of one singular value, 
with 



F R (X) 



dxp R (x 



o 



Actu ally , the singular values of a random matrix are not independent because of level repul 
sion 
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561 ] . It implies a zero probability for degenerate singular spaces: the singular values 
tend to keep away from each other. Thus, the probability density function pf (A) and the 
distribution function F R (X) of the first singular value Af cannot be deduced simply from 
p R (X). 

Once again, we will refer to RMT and particularly to the theoretical studies dealing with 



the first eigenvalue of the autocorrelation matrix K R K R t, i.e 

1 2 



R 



random matrix (multiple scattering regime) 



Af 



R 



4 + 



is given by 



Z + olNs 
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If K R is a "classical" 
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60|: 



(16) 



where Z is a random variable whose probability density function is a complicated law, 
known as the Tracy Widom distribution. It is an asymmetric bell curve, uncentered and of 
infinite support 
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]. No analytic expression is available, nevertheless it can be estimated 
numerically. In our case, we are rather interested in the probability density function pf of 
the first singular value Af . pf has been estimated by generating numerically "classical" 
random matrices and is displayed in Figfj^a), for A = 32 and A = 100. Although the 
quarter circle law pqc is of bounded support, pf is of infinite support. Indeed, pqc is only 
an asymptotic law, i.e valid for A — > oo. For a random matrix of finite dimension, Af 
has a non zero probability of exceeding A ma x = 2. Nevertheless, pf (A) narrows with A: 
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the variance of Af decreases when N grows, and its expected value has been shown to be 
bounded by A max = 2. 




Classic, N=32 

— Hankel, N=32 
■■■■ Classic, N=100 
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(a) (b) 

FIG. 9: (a) Probability density function pf of the first singular value Af estimated 
numerically. The multiple scattering case corresponds to continuous (N = 32) and 
dotted (N = 100) curves. The single-scattering case corresponds to dashed (N = 
32) and dash-dotted (N = 100) curves, (b) Distribution function F-[ i (\) of the 
first singular value Af . The multiple-scattering case still corresponds to continuous 
(N = 32) and dotted (N = 100) curves. The single-scattering case corresponds to 
dashed (N = 32) and dash-dotted (N = 100) curves. 



The relevant quantity for the detection issue is the distribution function of the first 
singular value Af . Ff' is the primitive of pf , 

F*(\) = p|Af<A} = J\xpf(x). 

In terms of target detection, the probability of false alarm is PFA(a) = 1 — F^(a). By 
setting an acceptable PFA, we can determine the corresponding threshold above which a 
target is said to be detected. The distribution functions F^, for N = 32 and N = 100, 
are plotted in Figj9](b). In the multiple scattering regime ("classical" random matrix), we 
observe that, for N = 32 and N = 100, Ff (A = 2) ~ 0, 99: if the PFA is fixed at 1%, a = 2 
is the corresponding threshold. 

This is no longer the case if the propagation operator behaves as a Hankel matrix (single 
scattering regime): the threshold is much higher, making the detection more difficult, for the 
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same PFA. Theoretical results indicate that the expected value of the first singular value 
grows as y^\og(N) 61]. This is illustrated in Figj9j pf and Ff 1 have been calculated by 
generating numerically random Hankel matrices for N = 32 and N = 100. As N is enlarged 
the curves are shifted towards right, which implies a larger probability of false alarm for a 
given detection threshold a. Other theoretical studies have shown that the fluctuations of 
Af are in 



bounded 



erior to y/\og N [62| and that, in the gaussian case, the variance of Af remains 



In a real situation, single and multiple scattering coexist; what is the relevant distribution 
function F^l If multiple scattering dominates (i.e., the target depth is larger than a few 
mean-free paths) or if the perturbation K R is additive uncorrelated noise, the probability 
of false alarm PFA(a) = 1 — F^(a) is deduced considering the distribution function F± 
obtained for a "classical" random matrix. If no a priori information is available regarding 
the scattering medium, we have to calculate Ff for a "classical" random matrix and for a 
random Hankel matrix, deduce two values for the probability of false alarm 1 — Ff(a), and 
keep the highest. 

In the view of applications (e.g., non-destructive evaluation, target detection etc.) an 
acceptable probability of false alarm P is set first, hence the detection threshold : 



a 



Ff^l-Po) (17) 



Once the usual normalization is performed, if Ai > a a target is detected at time T and 
frequency /. Conversely, if Ai < a one cannot conclude about the possible presence of a 
target. The application of this detection criterion has been performed experimentally in a 



recent study 64]. 



The detection criterion (Eq JTTl) can be used to estimate the performance of the DORT 
method in a random medium (or in the presence of noise). Let a T and a\ be the power of 
signals associated with the target and the random contribution. If the first singular value 
\i of K is associated to the target echo, then the expected value of Ai is Af . In Appendix 
|D| it is shown that Af = NaT, so 

E{Ax} = NaT ■ 

We also show in Appendix [D] that the quadratic mean of the singular values can be expressed 
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as 



N 

P =i 



Upon normalization, we have 



r 2 



JN ° T ~ —VN , for a 2 T « al 

a T + °R a R 



and var jAi j — \- Even though the echo of the target may be very weak compared to the 
scattering contribution, DORT can be expected to detect it if Ai > a, i.e., 

•* > f " (1 - p °>. (18) 
0R Vn 



The most favourable situation is when the quarter-circle law is valid (multiple scattering 
regime or additive white noise, and iV >> 1) : in that case a can be fixed to 2, which 
corresponds for N = 32 to a probability of false alarm P ~ 1%, and typically the target 
will be detected if 

°t 2 

In other words, the weakness of the target may be compensated by an increasing number 
of array elements: the performance of the DORT method is improved as the square-root of 
the number of independent channels. However in the single scattering regime, the expected 



value of Af grows as i/l°g N [6l(]. For a given Pq, we can assume that a increases also 
as y/logN. Hence, the performance of the D.O.R.T method grows as J j^y, which is 
significantly slower than V~N . 



VI. CONCLUSION 



In this article, the distribution of singular values of the array response matrix K in a 
random medium has been investigated, in a backscattering configuration. Once a judicious 
renormalization is achieved, the distribution p(A) obtained experimentally with ultrasonic 
waves in scattering samples is in very good agreement with theoretical predictions. Inter- 
estingly, p(A) is shown to differ significantly in the single and multiple scattering regimes. 
When multiple scattering dominates, as long as spatial correlations between matrix entries 
are negligible, p(X) is found to approach the quarter circle law. Correlations between ma- 
trix entries can also be taken into account: in that case, p(A) can be calculated following 
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the method proposed by Sengupta and Mitra[29|]. On the contrary, single scattering contri 



butions exhibit a different behaviour: whatever the realisation of disorder, a deterministic 
coherence persists along each antidiagonal of the matrix K. As a consequence, p(X) no longer 
follows the quarter- circle law and the singular spectrum of K becomes analogous to that of 
a Hankel matrix. These results have been applied to the detection of a target embedded 
in a random scattering media. Once the matrix is renormalized, knowing the distribution 
p(X) allows one to define a rigorous detection criterion based on the strongest singular value, 
which is expected to be associated to the target. The perspectives of this work are many. 
The results could be applied to all fields of wave physics where coherent transmit /receive 
arrays are available for imaging and detection (e.g., non destructive testing of scattering 
materials, underwater acoustics, landmine detection, seismology, radar/sonar, etc.) 
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Appendix A: 

The time-window length At is chosen so that signals associated with the same scattering 
path within the medium are in the same time window. For the sake of simplicity, here we only 
deal with single scattering paths but we would obtain the same results if multiple-scattering 
paths were considered. 

To calculate At, we have to find the two single-scattering paths, associated to the same 
scattering event, for which the difference of travel length is the highest. To that aim, the 
directivity of transducers has to be taken into account. The major part of the energy 
is transmitted (received) towards (from) scatterers located into a cone whose aperture is 
2@max- When the response is considered, only the scatterers contained inside the volume 
common to the directivity cones of elements % and j have to be considered. 

At each depth Z s , we have to optimize the lateral position X s of the scatterer which 
results in the largest difference of travel length between two single-scattering paths. These 
paths are associated with two source/receiver couples and (l,m). Thus, we have to 
optimize simultaneously the position X s , and the couples and (l,m). This issue can be 
simplified because we have to deal with single-scattering paths. Because of the equivalence 
of the forward and return waves, source and receiver of each couples are in fact identical. 
So, we have to find the elements % and / which are respectively the furthest and the nearest 
elements from the scatterer whose coordinates are (X S ,Z S ), X s remaining unknown. The 
result of this optimization differs according to the depth Z s of the scatterer : 
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FIG. 10: Choice of the appropriate time-window length At. The single scattered paths s and s' are 
associated with the reflection on the same scatterer. Two situations can occur. The case z = 
when < Zn m , corresponds to the dashed arrows. The two paths s\ and s'i have been chosen 

(2) (2) 

so that the difference of travel length is the largest. In the case z = Z s (with Z s > Zu m ), 
continuous arrows represent the two paths S2 and s' 2 for which the difference of travel length is the 
largest. 

• Z s < Z Um = y cos9 max : the two paths Si and s[ corresponding to the largest difference 
of travel length are shown with dashed arrows in FigHni si is the longest path and 
contributes to the signal ku. It is linked with a scatterer which is located on the top 
generatrix of the directivity cone of the first transducer. s[ is the shortest path. It is 
associated with the p th element of the array and thus contributes to the signal k pp . The 
transducer p is at the same transverse position X s as the scatterer. So, for Z s < Zn m , 
At is then given by 

(Al) 

• Z s > Zi im : the two paths S2 and s' 2 are depicted with continuous arrows in FigJTUl As 
previously, S2 is the longest path and contributes to the signal fen. It is still linked 
with a scatterer which is located on the top generatrix of the directivity cone of the 



At 



2Z„ 



COS 6 r , 
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first transducer. s' 2 is the shortest path. It is associated with the N th element of the 
array and thus contributes to the signal Unn- For Z s > Zu m , At is hence given by 

Z a 



At= 2 - 

c 



Z s tan 9 mnx . + D) 



cos max 

where D is the array aperture. When Z s — > oo, Eq JA2l becomes 

2D Sin 6 mnx 



Z? 



(A2) 



lim At 



(A3) 



In practice, we have to assess At by considering the maximum depth Z max = a + L. For 
the experiment described in SecJIIIj Z max = 65 mm. Considering an aperture angle 9 max = 
27.5 deg, we obtain a value At ~ 10/is. For the experiment described in SecJIVt Z max = 150 
mm. Thus the value At should be fixed to 14/is. Nevertheless, given the finite width of 
the agar gel sample (8 cm), there is no need to use such long temporal windows and a 
time-window length At = is also considered. 



Appendix B: 

As seen in Sec JIII Cj the coherent backscattering effect arises in the diagonal of matrix K: 
the variance of coefficients ka is twice that of off-diagonal elements. We want to estimate 
the influence of the coherent backscattering effect on the singular values distribution. To 
that aim, we will investigate its influence on the statistical properties of the autocorrelation 
matrix. 

Let us first consider a random matrix R of dimension N x N. We assume that the coef- 
ficients of the matrix R are complex gaussian random variables i.i.d, with a zero- mean and 
a variance of 1/N. The theoretical singular values distribution of a random matrix R (the 
quarter circle law) is deduced directly from the eigenvalues distribution of the autocorrela- 
tion matrix A = RR^ (the so-called Marcenko-Pastur law [l|), since singular values of R 
correspond to the square root of eigenvalues of A = RR^. We now focus on the statistical 
properties of A. The entries a/ m of matrix A are given by 

N 

aim = ^r lp r* mp (Bl) 
P =i 
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Let us calculate the mean and the variance of coefficients ai m . (ai m ) is given by 



N 

(aim) = ^2(ri P r* mp ) 



(aim) 



iil^m 

/I |2\ 



^2 p=1 (\ri p \ ) = 1 otherwise 
The mean of coefficients ai m is not nil only for diagonal elements 

(aim) = Sim (B2) 

{\a>im\ 2 ) can be developed as 

N N 

(\ a lm\ ) = ( r lp r mp r lq r mq) 
p=l q=l 

Using the moment theorem, (|a; m | 2 ) becomes 

N N N N 



,2 



) = Yl Yl ( ri P r *m P ) ( r lq r mq) (^h) ( r *mp r mq) (B3) 



p=l q=l p=l q=l 

v / s 



|(a im >| 2 var[a im ] 
To calculate var [ai m ] , we use the fact that 

We obtain 

1 N N 1 
var [a im ] = S ™ = n (B4) 

p=l q=l 

Finally, the off-diagonal entries of A are complex random variables with zero mean and 
variance 1/N. The diagonal elements of A are also complex random variables but with a 
mean equal to 1 and a variance of 1/N. 

Now that we have calculated the mean and variance of the autocorrelation coefficients a% m 
built from a "classical" random matrix R, we focus on the effect of coherent backscattering. 
To that aim, we consider the matrix K which is obtained experimentally in the multiple 
scattering regime. This matrix is built from the array response matrix K: 

fey = ^ = (B5) 

N Z^p=l Z^iq=l \ fh pq\ 
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Eq JB5l is another expression for the renormalization presented in EqJU By construction, we 
have 

vEEN =i (b6) 



t=i j=i 

Note that the random matrix R verifies the same property. But, as pointed out in Sec JIII C| 
the variance of elements of K is not constant. In particular, the matrix K exhibits a doubled 
variance along its diagonal because of the coherent backscattering phenomenon. To estimate 
its effect on the singular value spectrum, we model K as a matrix whose coefficients k\j are 
complex gaussian random variables independently distributed with mean zero and a variance 
of,- defined as 




Eq JB6l leads to 



1 N N 1 

i\T2 °V<1 = ft 

p=l q=l 

p=l 

N(N + 1) 2 J_ 

iV 2 a TV' 



which yields 



of- is finally given by 



a 2 = — ^- (B8) 
TV + 1 V ; 



a% = i±% (B9) 



The variance of elements fc^- is equal to j^tj for off-diagonal elements (i ^ j) and twice for 
diagonal elements (z = j). 

Let us focus now on the autocorrelation matrix B = KKt Its coefficients are defined as 

N 

km = Yl KK (BIO) 

p=l 
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The mean of coefficients km is given by 



iV 



(km) = ^2(klpK 
p=l 



mp 



(km) 



0, if I y£ m 



Hp 



p^l u lp N+l 



jj—j = 1, otherwise 



We obtain the same mean for km as for a\ m (Eq JB2l) 

(km) = 8lm 

As to the variance of coefficient k m , it can be expressed as 

N N 



(Bll) 



var [km] 

var [bi m } 

var [6 Im ] 

var [k m ] 
var [fy m ] 



p=l q=l 



mq 



t 2 r) 

7 mq°Pq 



N 



E 

9=1 



2 2 



(iV 



1 * 



1 



(A+l) 2 
A + 2 + 5 lm 



(A + 2 + 5 lm ) 



(B12) 



(A+l) 2 

The statistical properties of the autocorrelation matrices A and B are summed up in TabJT] 
The coherent backscattering effect increases the variance of the autocorrelation coefficients, 



"Classical" random matrix R 
A = RRt 


Experimental matrix K 
B = KKt 


(aim) = km 


(km) = Sim 


var [a im ] = ± 


var [Oim\ — (tv+i)2 



TABLE I: Statistical properties of autocorrelation matrices A and B, depending on whether the 
initial matrix exhibits the coherent backscattering effect or not. 



compared to the classical case. Thus, the eigenvalues spectrum of the autocorrelation matrix 
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and also the singular values distribution of K should be modified by the coherent backscat- 
tering effect. But, for N >> 1, the difference between r^fp and t^^s becomes 
negligible. For instance, the relative error is about 3% for N = 32. So, the doubled variance 
of diagonal elements of K does not perturb significantly the singular values distribution. 



Appendix C: 

We investigate the influence of reciprocity on the statistical properties of the autocorre- 
lation matrix. We will consider first the case of a "classical" random matrix and study the 
statistical properties of its autocorrelation matrix. The case of a symetric random matrix 
will be studied afterwards. The comparison of the results obtained in both cases will allow 
us to explain and quantify the influence of symmetry on the singular values distribution. 
As we will see, the symmetry implies correlations between the diagonal elements of the 
autocorrelation matrix. 

Let us first consider a random matrix R of dimension N x N. We assume that the 
coefficients of the matrix R are complex gaussian random variables i.i.d, zero-mean and 
with variance 1/N. The entries a\ m of the autocorrelation matrix A = RR^ are given by 



iV 



0>lr, 



} j r lp r mp 
p=l 



Let us calculate the correlation coefficient 9^ between diagonal elements an and a mm , which 
is defined as 

A _ (a u a* mm ) - \(a u )\ 
The correlation term (aua* mm ) can be developed as 

N M 

( a U a mm) = { ri P r lp rm i r mq) 

Using the moment theorem, the last equation becomes 

N M N M 

(a a a* mm ) = Yl ( ri P r t) { r W*mq) + ( ri P r *mq) (^^p) (C2) 

p=l q=l p=l q=l 

V v ' 

The second sum can be calculated using the fact that 

(n P r* mq ) = %^ (C3) 
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which yields 

N M 



{a u a* mm } - \ (a a }\ 2 = 



3lm&pq 

p=l q=l 



N 2 



(a u a* mm ) - \ (a u }\ 2 = 

Using the fact that var [an] = (EqjB4]), the correlation coefficient 9fa„(Eq JCll is given by 

®L = Sim (C4) 

This last equation means that the diagonal elements an of A are totally decorrelated from 
each other. 

We now consider the case of a random but symmetric matrix like K. Its coefficients 
are complex gaussian random variables identically distributed, with mean zero and variance 
1/N. The only difference with a "classical" random matrix R is the fact that ky = kji. 
We neglect here the coherent backscattering effect (see Appendix [B]) in order to focus on the 
influence of reciprocity. 

Let us study the statistical properties of the autocorrelation matrix B = KKt One can 
show that symmetry has no effect on the mean and the variance of b\ m 

(km) = (aim) = -j^ 

var [bim] = var [a im ] = — 

As done previously for the matrix A, we define the correlation coefficient Qf m between 
diagonal elements bu and b mm 

_ ( b ub* mm ) - | (Ml 2 
var [bu\ 

The correlation term (bab mm ) can be developed as 

N M 



^U^mm) ~ \klpklpkmqk m 
P=l 9=1 



Using the moment theorem, the last equation becomes 

N M N M 



ihlKnm) ~ (^P lp/ \kmqk mq \ + y^lp^mq/ \k mq k lp 



p=l q=l p=l q=l 

v ' 

\(hl)\ 2 
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But, because ki p = k p i, we have 

{hK) = ^ + 6 -^l (C6) 

We obtain: 

( b ub* mm ) - \ (bu)\ 2 = -j^- + 
Using the fact that var [bu] = jj, the correlation coefficient 6^(Eq lC5l) is finally given by 

©f m = Slm + Jf (C7) 

This last equation means that a correlation exists between the diagonal elements bu of B, 
due to reciprocity. The correlations between diagonal elements of matrices A and B are 
compared in TabHTl Even if this correlation is not too large (1/iV ~ 3% with N = 32), it 



"Classical" random matrix R 
A = RRt 


Symmetric random matrix K 
B = KK7 




®hn = km + n 



TABLE II: Correlation between diagonal elements of autocorrelation matrices A and B, depending 
on whether the initial matrix is symmetric or not. 



has an influence on the eigenvalue spectrum of the autocorrelation matrix, hence on the 
singular values distribution of the matrix K. This explains partly the deviation from the 
quarter circle law, which we pointed out numerically (see Figj6]^b)). Note that for iV >> 1, 
the influence of symmetry should vanish. 



Appendix D: 

We assume that the matrix K T associated to the target is of rank 1. Af is the only 
non-zero singular value of K T . The trace of the autocorrelation matrix K T K T ^ can then be 
expressed as 

Trace 



N 

K T K Tt ]=£(Aj) 2 =(Af) 2 (D1) 

p=i 

a\ is the power of signals associated with the target : 



j N N 
i=X j=l 
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The trace of K T K T ^ can also be expressed as 

Trace 



K T K Tt 



N N 



Trace 

From Eq JDll and Eq JD2l we deduce 



K T K Tt 



EE 1*5 

i=i j=i 

N 2 a 2 T (D2) 



\\ = Na T (D3) 

Now, we consider the case of the matrix K. This matrix is the sum of two contributions 
(EqlE]) : 

• K T which corresponds to the direct echo reflected by the target. 

• K R which corresponds to the response of the random medium (or to additive noise). 
The trace of the autocorrelation matrix KK^ is given by 

N 

Trace [KK f ] = K ( D4 ) 
P =i 

We note a\ the mean power of signals linked with the random contribution: 

. N N 
i=l j=l 

The trace of KK* can be also expressed as 



N N 
i=l j=l 




N N 

EE(*5h 

i=l j=l 




N N 

eei^i 2 


N N N N 

+EEi*sr+ 2 EE*[^r] 


i=l j=l 

V v ' 

iV 2 c4 


i=l j=l i=l i=l 

" v ' 



Assuming that Nar » ctr and Nctr » err (assumptions verified a posteriori with Eq JT9|) . 
we can neglect the third sum because (JR, [kjjk^j*] ) = and std Yli=i ^2f=i ^ [kfj^ij*] = 
NcttCr « N 2 o~t, N 2 a\. And we finally obtain 

Trace [KK f ] = N 2 {o 2 T + a\) (D5) 
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From Eq lD4l and EqjD5j we can deduce an expression for the quadratic mean of the singular 
values 

^E A p= \fn^+^) ( D6 ) 
P =i 



41 



